

 ########################################################################
 ## INCENTIVIZING ANTICORRUPTION REFORM  - APPENDIX REPLICATION FILES  ##
 ########################################################################

 ## Last modified: NOV, 2024

 # Clean
 rm(list=ls())

 # Set directory
 setwd("")

 # Packages
 library(ggplot2)
 library(dplyr)
 library(ggpubr)
 library(did)
 library(doBy)
 library(stargazer)
 library(MASS)
 library(pscl)
 library(boot)
 library(reshape2)
 library(stringr)
 library(marginaleffects)
 library(readxl)


 ## CONTENT:
 # 1) Summary statistics (appendix section 8)
 # 2) Descriptive plots of ACIs (appendix section 9)
 # 3) Multiple periods difference-in-differences (appendix section 10)
 # 4) Additional regression tables (appendix section 11)
 # 5) Full legislative term models (appendix section 12)
 # 6) Party analysis (appendix section 13)
 
 
 
 ################################################ 1) SUMMARY STATISTICS

 ## TABLE A.6 - Summary statistics (2015-2017 models)

 # Import
 leg_did <- readRDS(file = "leg_DID.rds")
 
 # Choose vars
 vars <- leg_did[, c("acl_total", "sna_related", "reelection", "pres_coal",  
                     "gov_coal", "margin_victory", "pop_log", "per_rural", 
                     "corr_perception", "l_dep_budget")]
 
 # Table
 stargazer(vars, summary.stat = c("n", "mean", "sd", "min", "max"))
 
 
 ## TABLE A.7 - Summary statistics (full legislative term models)
 
 # Import
 legs <- readRDS(file = "full_term.rds")

 # Create numeric variables
 legs$smd <- 0
 legs$smd[legs$district_type=="SMD"] <- 1
 legs$pr <- 0
 legs$pr[legs$district_type=="PR"] <- 1
 
 # Choose vars
 vars <- legs[, c("acl", "reelection",  "smd", "pr", "gov_coal")]
 
 # Table
 stargazer(vars, summary.stat = c("n", "mean", "sd", "min", "max"))
 
 
 
 
 
 
 ################################################ 2) DESCRIPTIVE PLOTS FOR ACIs
 
 
 #### FIGURE A.3  - Anticorruption initiatives (ACIs) in time
 
 # Import
 plot1 <- read.csv("Figure A3.csv")

 ## PLOT
 g1 <- ggplot(plot1, aes(x=year, y=ACI)) + 
   geom_line(color="darkgray", size=0.5) + 
   theme_light() +
   ggtitle("Anticorruption initiatives (ACIs) in time") +
   xlab("") + 
   ylab("ACIs") +
   ylim(0, 200) +
   scale_x_continuous(breaks=c(2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020)) +
   theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))  +  
   theme(plot.margin=unit(c(.5, .5 , .8, .5),"cm"),
         plot.title = element_text(hjust = 0.5),
         panel.grid.major = element_blank(),
         panel.grid.minor = element_blank(),
         axis.title = element_text(margin = margin(20, 20, 20, 20)),
         axis.title.x.bottom = element_text(margin = margin(0, 0, 1, 0)),
         axis.title.y.left = element_text(margin = margin(0, 15, 0, 1))) +
   geom_vline(xintercept = 2014.2, linetype="dotted", color = "black", size=0.5) +
   annotate(geom = "text", x = 2014, y = 150, label = "White House Scandal", color = "black", angle = 90, size = 3) +
   geom_vline(xintercept = 2014.9, linetype="dotted", color = "black", size=0.5) +
   annotate(geom = "text", x = 2014.7, y = 150, label = "Oceanografia Scandal", color = "black", angle = 90, size = 3) +
   geom_vline(xintercept = 2015.3, linetype="dotted", color = "black", size=0.5) +
   annotate(geom = "text", x = 2015.5, y = 156, label = "SNA constitutional reform", color = "black", angle = 90, size = 3) +
   geom_vline(xintercept = 2016.5, linetype="dotted", color = "black", size=0.5) +
   annotate(geom = "text", x = 2016.3, y = 155, label = "SNA secondary reforms", color = "black", angle = 90, size = 3)
 g1
 
 
 
 
 #### FIGURE A.4  - Anticorruption initiatives (ACIs) and term limit status

 # Import 
 plot2 <- read.csv("Figure A4.csv")
 
 # Categorical variables
 plot2$reelection[plot2$reelection==0] <- "Term-limited"
 plot2$reelection[plot2$reelection==1] <- "Reelection incentives"
 
 # Order factor levels
 plot2$state <- factor(plot2$state, levels = c("Chiapas", "Chihuahua", "Morelos", "Queretaro", 
                                             "Quintana Roo", "Campeche", "Colima", "Sinaloa",
                                             "Tlaxcala", "Baja California Sur", "Guanajuato", 
                                             "Guerrero", "Hidalgo", "Jalisco", "Nuevo Leon",
                                             "Tabasco", "Zacatecas", "Aguascalientes",
                                             "Durango", "Nayarit", "Michoacan", "Oaxaca",
                                             "CDMX", "Puebla", "Tamaulipas", "Baja California", 
                                             "Sonora", "Coahuila")) 
 
 
 # Exclude Mexico, SLP, Yucatan, and and Veracruz due to incomplete data
 plot2 <- plot2[-c(29,30,47,48,59,60,61,62) , ]
 
 ## PLOT
 g2 <- ggplot(plot2, aes(fill=reelection, y=ACL, x=forcats::fct_infreq(state))) + 
   geom_bar(position="stack", stat="identity") +
   theme_light() +
   ggtitle("Anticorruption initiatives and term limit status") +
   xlab("") + 
   ylab("ACIs") + 
   coord_flip() + 
   theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))  +  
   theme(legend.title=element_blank(),
         #legend.position = c(0.85, .925),
         legend.position = c(0.8, .1),
         legend.background = element_rect(fill = "white", color = "gray"),
         panel.grid.major = element_blank(), 
         #panel.grid.minor = element_blank(),
         legend.key.size = unit(0.3, 'cm'),
         plot.title = element_text(hjust = 0.5),
         plot.margin=unit(c(0.5, 1 , .8, .5),"cm"),
         axis.title.x.bottom = element_text(margin = margin(5, 0, 1, 0)),
         axis.title.y.left = element_text(margin = margin(0, 5, 0, 1)))
 g2 <-  g2 + scale_fill_manual(values = c("#696969", "#BEBEBE"))
 g2
 
 
 
 
 #### FIGURE A.5  - Subject matter of anticorruption initiatives
 
 topics <- readRDS(file = "ACI_content.rds")
 
 # Create count
 topics$total <- 1
 
 cats <- topics[, c("total", "topics")]
 
 # Note and re-order
 table(cats$topics, exclude = NULL)
 
 cats$topics <- trimws(cats$topics)
 cats <- na.omit(cats)
 
 cats$topics <- ordered(cats$topics , levels = c(
   "Financial disclosures", 
   "Whistleblower protections", 
   "Electoral malfeasance",
   "Anticorruption institutions",
   "Appointments and background checks",
   "Misappropriation and financial corruption",
   "Trading in influence", 
   "Public procurement",
   "Awareness",
   "Normativity",
   "Sanctions for corruption",
   "Transparency and open government",
   "Audits and investigations",
   "SNA-related"))
 
 ## PLOT
 ggplot(cats, aes(y=total, x=topics)) + 
   geom_bar(position="stack", stat="identity", width = 0.9) +
   theme_light() +
   ggtitle("") +
   xlab("") + 
   ylab("") + 
   coord_flip() + 
   theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))  +  
   theme(legend.title=element_blank(),
         #legend.position = c(0.85, .925),
         legend.position = c(0.8, .1),
         legend.background = element_rect(fill = "white", color = "gray"),
         panel.grid.major = element_blank(), 
         #panel.grid.minor = element_blank(),
         legend.key.size = unit(0.3, 'cm'),
         plot.title = element_text(hjust = 0.5),
         plot.margin=unit(c(0.5, 1 , .8, .5),"cm"),
         axis.title.x.bottom = element_text(margin = margin(5, 0, 1, 0)),
         axis.title.y.left = element_text(margin = margin(0, 5, 0, 1)))
 
 
 
 
 
 
 ################################################ 3) MULTIPLE PERIODS DIFFERENCE-IN-DIFFERENCES
 

 ## TABLE A.8 - Estimates of the effect of reelection on SNA bill sponsorship (never-treated control)
 
 # Never treated
 nt <- att_gt(yname = "sna_related",
              tname = "time_var",
              gname = "treat_g",
              idname = "legislator_ID",
              xformla = NULL,
              allow_unbalanced_panel = TRUE,
              data = leg_did,
              est_method = "reg", 
              panel = FALSE,
              clustervars = c("treat_clustvar"), 
              bstrap = FALSE,
              control_group = c("nevertreated"))
 
 
 ## Simple aggregation
 agg.simple <- aggte(nt, type = "simple")
 summary(agg.simple)
 
 ## Specific groups
 agg.gs <- aggte(nt, type = "group")
 summary(agg.gs)
 ggdid(agg.gs) + theme_light() 
 


 ## FIGURE A.6 - Group-time average treatment effects (not-yet-treated)
 nyt <- att_gt(yname = "sna_related",
               tname = "time_var",
               gname = "treat_g",
               idname = "legislator_ID",
               xformla = NULL,
               allow_unbalanced_panel = TRUE,
               data = leg_did,
               est_method = "reg", 
               panel = FALSE,
               clustervars = c("treat_clustvar"), 
               bstrap = FALSE,
               control_group = c("notyettreated"))
 
 ggdid(nyt, theming  = FALSE, legend = FALSE) +
   theme_light() +
   ggtitle("") + 
   ylab("SNA related bills") +
   #ylim(-0.25, 1.5) + 
   theme(plot.margin=unit(c(1, 1 , .8, .5),"cm"),
         panel.grid.minor = element_line(size = 0.1),
         panel.grid.major = element_line(size = 0.1),
         axis.title = element_text(margin = margin(20, 20, 20, 20)),
         axis.title.x.bottom = element_text(margin = margin(10, 0, 1, 0)),
         axis.title.y.left = element_text(margin = margin(1, 15, 0, 1))) +
   scale_fill_discrete(labels=c('Term-limited', 'Reelection incentives'))
 
 
 
 ## FIGURE A.7 - Group-time average treatment effects (never-treated)
 ggdid(nt, theming  = FALSE, legend = FALSE) +
   theme_light() +
   ggtitle("") + 
   ylab("SNA related bills") +
   #ylim(-0.25, 1.5) + 
   theme(plot.margin=unit(c(1, 1 , .8, .5),"cm"),
         panel.grid.minor = element_line(size = 0.1),
         panel.grid.major = element_line(size = 0.1),
         axis.title = element_text(margin = margin(20, 20, 20, 20)),
         axis.title.x.bottom = element_text(margin = margin(10, 0, 1, 0)),
         axis.title.y.left = element_text(margin = margin(1, 15, 0, 1))) +
   scale_fill_discrete(labels=c('Term-limited', 'Reelection incentives'))
 
 
 
 
 
 ################################################ 4) ADDITIONAL REGRESSION TABLES
 
 ## TABLE A.9 - Reelection incentives, district type, and membership to the governor’s coalition (binary outcomes)

 # SNA DUMMY
 leg_did$sna_d <- leg_did$sna_related
 leg_did$sna_d[leg_did$sna_d>0] <- 1

 # ACI DUMMY
 leg_did$aci_d <- leg_did$acl_total
 leg_did$aci_d[leg_did$aci_d>0] <- 1
 
 # Change to numeric
 leg_did$smd <- NA
 leg_did$smd[leg_did$district_type=="SMD"] <- 1
 leg_did$smd[leg_did$district_type=="PR"] <- 0

 
 # 1) SNA state and period FE (no controls)
 mod1 <- lm(sna_d ~ reelection + factor(state) + factor(year), data = leg_did)
 
 # State CSE
 summary(mod1, cluster = c("state"))
 mod_1_se <- summary(mod1, cluster = c("state"))$coefficients[,"Std. Error"]
 
 
 # 2) ACI state and period FE (no controls)
 mod2 <- lm(aci_d ~ reelection + factor(state) + factor(year), data = leg_did)
 
 # State CSE
 summary(mod2, cluster = c("state"))
 mod_2_se <- summary(mod2, cluster = c("state"))$coefficients[,"Std. Error"]
 
 
 # 3) SNA state and period FE (controls)
 mod3 <- lm(sna_d ~ reelection + gov_coal + pres_coal + margin_victory + 
              pop_log + per_rural + corr_perception +  l_dep_budget + 
              factor(state) + factor(year), data = leg_did)
 
 # State CSE
 summary(mod3, cluster = c("state"))
 mod_3_se <- summary(mod3, cluster = c("state"))$coefficients[,"Std. Error"]
 
 
 # 4) ACI state and period FE (controls)
 mod4 <- lm(aci_d ~ reelection + gov_coal + pres_coal + margin_victory + 
              pop_log + per_rural + corr_perception + l_dep_budget + 
              factor(state) + factor(year), data = leg_did)
 
 # State CSE
 summary(mod4, cluster = c("state"))
 mod_4_se <- summary(mod4, cluster = c("state"))$coefficients[,"Std. Error"]
 
 
 # 5) SNA INTERACTION REELEC X DISTRICT TYPE
 mod5 <- lm(sna_d ~ reelection*smd + gov_coal + pres_coal + margin_victory + 
              pop_log + per_rural + corr_perception +  l_dep_budget + 
              factor(state) + factor(year), data = leg_did)
 
 # State CSE
 summary(mod5, cluster = c("state"))
 mod_5_se <- summary(mod5, cluster = c("state"))$coefficients[,"Std. Error"]
 
 
 # 6) ACI INTERACTION REELEC X DISTRICT TYPE
 mod6 <- lm(aci_d ~ reelection*smd + gov_coal + pres_coal + margin_victory + 
              pop_log + per_rural + corr_perception +  l_dep_budget + 
              factor(state) + factor(year), data = leg_did)
 
 # State CSE
 summary(mod6, cluster = c("state"))
 mod_6_se <- summary(mod6, cluster = c("state"))$coefficients[,"Std. Error"]
 
 
 # 7) SNA INTERACTION REELEC X GOVS
 mod7 <- lm(sna_d ~ reelection*gov_coal + pres_coal + margin_victory + 
              pop_log + per_rural + corr_perception +  l_dep_budget + 
              factor(state) + factor(year), data = leg_did)
 
 # State CSE
 summary(mod7, cluster = c("state"))
 mod_7_se <- summary(mod7, cluster = c("state"))$coefficients[,"Std. Error"]
 
 
 # 8) ACI INTERACTION REELEC X GOVS
 mod8 <- lm(aci_d ~ reelection*gov_coal + pres_coal + margin_victory + 
              pop_log + per_rural + corr_perception +  l_dep_budget + 
              factor(state) + factor(year), data = leg_did)
 
 # State CSE
 summary(mod8, cluster = c("state"))
 mod_8_se <- summary(mod8, cluster = c("state"))$coefficients[,"Std. Error"]
 
 
 
 ## Table
 stargazer(mod1, mod2, mod3, mod4, mod5, mod6, mod7 , mod8 ,
           title="Reelection and anticorruption bill sponsorship",
           se = list(mod_1_se, mod_2_se, mod_3_se, mod_4_se, 
                     mod_5_se, mod_6_se, mod_7_se, mod_8_se),
           align=TRUE, digits = 2, 
           digits.extra = 0, 
           dep.var.labels=c("ACI", "SNA"),
           covariate.labels=c("Reelection incentives"))
 

 
 ## TABLE A.10 - Reelection incentives, district type, and membership to the governor’s coalition (negative binomials)

 
 # 1) SNA state and period FE (no controls)
 mod1 <- glm.nb(sna_related ~ reelection + factor(state) + factor(year), data = leg_did)
 
 # State CSE
 summary(mod1, cluster = c("state"))
 mod_1_se <- summary(mod1, cluster = c("state"))$coefficients[,"Std. Error"]
 
 
 # 2) ACI state and period FE (no controls)
 mod2 <- glm.nb(acl_total ~ reelection + factor(state) + factor(year), data = leg_did)
 
 # State CSE
 summary(mod2, cluster = c("state"))
 mod_2_se <- summary(mod2, cluster = c("state"))$coefficients[,"Std. Error"]
 
 
 # 3) SNA state and period FE (controls)
 mod3 <- glm.nb(sna_related ~ reelection + gov_coal + pres_coal + margin_victory + 
                  pop_log + per_rural + corr_perception +  l_dep_budget + 
                  factor(state) + factor(year), data = leg_did)
 # State CSE
 summary(mod3, cluster = c("state"))
 mod_3_se <- summary(mod3, cluster = c("state"))$coefficients[,"Std. Error"]
 
 
 # 4) ACI state and period FE (controls)
 mod4 <- glm.nb(acl_total ~ reelection + gov_coal + pres_coal + margin_victory + 
                  pop_log + per_rural + corr_perception + l_dep_budget + 
                  factor(state) + factor(year), data = leg_did)
 # State CSE
 summary(mod4, cluster = c("state"))
 mod_4_se <- summary(mod4, cluster = c("state"))$coefficients[,"Std. Error"]
 
 
 # 5) SNA INTERACTION REELEC X DISTRICT TYPE
 mod5 <- glm.nb(sna_related ~ reelection*smd + gov_coal + pres_coal + margin_victory + 
                  pop_log + per_rural + corr_perception +  l_dep_budget + 
                  factor(state) + factor(year), data = leg_did)
 
 # State CSE
 summary(mod5, cluster = c("state"))
 mod_5_se <- summary(mod5, cluster = c("state"))$coefficients[,"Std. Error"]
 
 
 # 6) ACI INTERACTION REELEC X DISTRICT TYPE
 mod6 <- glm.nb(acl_total ~ reelection*smd + gov_coal + pres_coal + margin_victory + 
                  pop_log + per_rural + corr_perception +  l_dep_budget + 
                  factor(state) + factor(year), data = leg_did)
 # State CSE
 summary(mod6, cluster = c("state"))
 mod_6_se <- summary(mod6, cluster = c("state"))$coefficients[,"Std. Error"]
 
 
 # 7) SNA INTERACTION REELEC X GOVS
 mod7 <- glm.nb(sna_related ~ reelection*gov_coal + pres_coal + margin_victory + 
                  pop_log + per_rural + corr_perception +  l_dep_budget + 
                  factor(state) + factor(year), data = leg_did)
 # State CSE
 summary(mod7, cluster = c("state"))
 mod_7_se <- summary(mod7, cluster = c("state"))$coefficients[,"Std. Error"]
 
 
 # 8) ACI INTERACTION REELEC X GOVS
 mod8 <- glm.nb(acl_total ~ reelection*gov_coal + pres_coal + margin_victory + 
                  pop_log + per_rural + corr_perception +  l_dep_budget + 
                  factor(state) + factor(year), data = leg_did)
 # State CSE
 summary(mod8, cluster = c("state"))
 mod_8_se <- summary(mod8, cluster = c("state"))$coefficients[,"Std. Error"]
 
 
 
 ## Table
 stargazer(mod1, mod2, mod3, mod4, mod5, mod6, mod7 , mod8 ,
           title="Reelection and anticorruption bill sponsorship",
           se = list(mod_1_se, mod_2_se, mod_3_se, mod_4_se, 
                     mod_5_se, mod_6_se, mod_7_se, mod_8_se),
           align=TRUE, digits = 2, 
           digits.extra = 0, 
           dep.var.labels=c("ACI", "SNA"),
           covariate.labels=c("Reelection incentives"))
 
 
 
 
 
 
 
 ################################################ 5) FULL LEGISLATIVE TERM MODELS
 
 ## TABLE A.11 - Reelection incentives, district type, and governor’s coalition (full legislative terms)
 
 # Create SMD and PR subsets
 smds <- subset(legs, legs$district_type=="SMD")
 pr <- subset(legs, legs$district_type=="PR")
 
 
 ## All (no FE)
 mod1 <- lm(acl ~ reelection, data = legs)
 summary(mod1)
 
 summary(mod1, cluster = c("state"))
 mod_1_se <- summary(mod1, cluster = c("state"))$coefficients[,"Std. Error"]
 
 # All (state FE)
 mod2 <- lm(acl ~ reelection + factor(state), data = legs)
 summary(mod2)
 
 summary(mod2, cluster = c("state"))
 mod_2_se <- summary(mod2, cluster = c("state"))$coefficients[,"Std. Error"]
 
 ## SMDs (no FE)
 mod3 <- lm(acl ~ reelection, data = smds)
 summary(mod3)
 
 summary(mod3, cluster = c("state"))
 mod_3_se <- summary(mod3, cluster = c("state"))$coefficients[,"Std. Error"]
 
 # SMDs (state FE
 mod4 <- lm(acl ~ reelection + factor(state), data = smds)
 summary(mod4)
 
 summary(mod4, cluster = c("state"))
 mod_4_se <- summary(mod4, cluster = c("state"))$coefficients[,"Std. Error"]
 

 ## PR (no FE)
 mod5 <- lm(acl ~ reelection, data = pr)
 summary(mod5)
 
 summary(mod5, cluster = c("state"))
 mod_5_se <- summary(mod5, cluster = c("state"))$coefficients[,"Std. Error"]
 
 # PR (state FE)
 mod6 <- lm(acl ~ reelection + factor(state), data = pr)
 summary(mod6)
 
 summary(mod6, cluster = c("state"))
 mod_6_se <- summary(mod6, cluster = c("state"))$coefficients[,"Std. Error"]
 
 
 ## INTERACTION SMD vs PR
 legs$smd <- NA
 legs$smd[legs$district_type=="SMD"] <- 1
 legs$smd[legs$district_type=="PR"] <- 0

 mod7 <- lm(acl ~ reelection*smd + factor(state), data = legs)
 summary(mod7)
 
 summary(mod7, cluster = c("state"))
 mod_7_se <- summary(mod7, cluster = c("state"))$coefficients[,"Std. Error"]
 
 
 ## INTERACTION GOVS
 mod8 <- lm(acl ~ reelection*gov_coal + factor(state), data = legs)
 summary(mod8)
 
 summary(mod8, cluster = c("state"))
 mod_8_se <- summary(mod8, cluster = c("state"))$coefficients[,"Std. Error"]
 
 
 
 stargazer(mod1, mod2, mod3, mod4, mod5, mod6, mod7, mod8,
           title="Reelection incentives and sponsored ACIs",
           se = list(mod_1_se, mod_2_se, mod_3_se, mod_4_se, mod_5_se, 
                     mod_6_se, mod_7_se, mod_8_se),
           align=TRUE, digits = 2, 
           digits.extra = 0, 
           dep.var.labels=c("ACIs"),
           covariate.labels=c("Reelection incentives"),
           column.labels   = c("Full sample", "SMD", "PR"),
           column.separate = c(2, 2, 2))
 
 
 
 
 
 ################################################ 6) PARTY ANALYSIS
 
 ## FIGURE A.8 - Anticorruption sponsorship by party and co-sponsorship status

 # Clean party variable 
 leg_did$party_clean <- leg_did$party
 leg_did$party_clean[leg_did$party_clean=="Alternativa Veracruzana"] <- "Other"
 leg_did$party_clean[leg_did$party_clean=="Chiapas Unido"] <- "Other"
 leg_did$party_clean[leg_did$party_clean=="Compromiso por Puebla"] <- "Other"
 leg_did$party_clean[leg_did$party_clean=="Mover a Chiapas"] <- "Other"
 leg_did$party_clean[leg_did$party_clean=="Partido de Baja California"] <- "Other"
 leg_did$party_clean[leg_did$party_clean=="Partido Conciencia Popular"] <- "Other"
 leg_did$party_clean[leg_did$party_clean=="Partido Alianza Ciudadana"] <- "Other"
 leg_did$party_clean[leg_did$party_clean=="Partido Duranguense"] <- "Other"
 leg_did$party_clean[leg_did$party_clean=="Partido Renovacion Sudcaliforniana"] <- "Other"
 leg_did$party_clean[leg_did$party_clean=="Partido Sinaloense"] <- "Other"
 leg_did$party_clean[leg_did$party_clean=="Partido Unidad Popular"] <- "Other"
 leg_did$party_clean[leg_did$party_clean=="Transformemos"] <- "Other"
 leg_did$party_clean[leg_did$party_clean=="UDC"] <- "Other"
 leg_did$party_clean[leg_did$party_clean=="PSD"] <- "Other"
 leg_did$party_clean[leg_did$party_clean=="Partido Socialdemocrata de Oaxaca"] <- "Other"
 leg_did$party_clean[leg_did$party_clean=="Partido Socialdemocrata Independiente"] <- "Other"
 leg_did$party_clean[leg_did$party_clean=="Partido Socialista "] <- "Other"
 leg_did$party_clean[leg_did$party_clean=="Primero Coahuila"] <- "Other"
 leg_did$party_clean[leg_did$party_clean=="PH"] <- "Other"
 leg_did$party_clean[leg_did$party_clean=="CONV"] <- "Other"
 leg_did$party_clean[leg_did$party_clean=="Partido Socialista"] <- "Other"
 leg_did$party_clean[leg_did$party_clean=="PRI-PVEM"] <- "PRI"
 
 table(leg_did$party_clean)
 
 # Omit NAs in party
 leg_did <- subset(leg_did , leg_did$party_clean !="")
 
 
 ### PLOT PARTY DESCRIPTIVES
 
 plot <- leg_did[, c("acl_total", "solo_author", "all_party", "party_clean")]
 
 plot <- summaryBy(acl_total + solo_author + all_party ~ party_clean, 
                   FUN=c(sum), data=plot, na.rm=TRUE, keep.names = TRUE)
 
 plot$party_clean <- factor(plot$party_clean, levels = c("PAN",
                                                         "PRI", 
                                                         "PRD", 
                                                         "MORENA",
                                                         "MC",
                                                         "PT",
                                                         "PVEM",
                                                         "PANAL",
                                                         "PES",
                                                         "Other"))
 plot <- na.omit(plot)
 
 p1 <- plot[, c("party_clean", "acl_total")]
 p1$group <- "ACIs"
 colnames(p1)[2] <- "total"
 
 p2 <- plot[, c("party_clean", "all_party")]
 p2$group <- "All party"
 colnames(p2)[2] <- "total"
 
 p3 <- plot[, c("party_clean", "solo_author")]
 p3$group <- "Single sponsor"
 colnames(p3)[2] <- "total"
 
 plot <- rbind(p1, p2, p3)
 
 
 ## Figure A.8
 ggplot(plot, aes(y=total, x=party_clean, fill = group)) + 
   geom_bar(position="dodge", stat="identity") +
   theme_light() +
   ggtitle("") +
   scale_fill_grey(start = .8, end = .2) +
   xlab("") + 
   ylab("ACIs sponsored by deputy") + 
   theme(legend.title=element_blank(),
         #legend.position = c(0.85, .925),
         legend.position = c(0.8, 0.8),
         legend.background = element_rect(fill = "white", color = "gray"),
         panel.grid.major = element_blank(), 
         #panel.grid.minor = element_blank(),
         legend.key.size = unit(0.3, 'cm'),
         plot.title = element_text(hjust = 0.5),
         axis.title.x.bottom = element_text(margin = margin(5, 0, 1, 0)),
         axis.title.y.left = element_text(margin = margin(0, 10, 1, 1))) 
 
 
 
 
 
 ## FIGURE A.9 Estimates of reelection incentives and anticorruption sponsorship by party 
 
 ## SNA
 pan <- subset(leg_did, leg_did$party_clean=="PAN")
 pri <- subset(leg_did, leg_did$party_clean=="PRI")
 mor <- subset(leg_did, leg_did$party_clean=="MORENA")
 prd <- subset(leg_did, leg_did$party_clean=="PRD")
 mc <- subset(leg_did, leg_did$party_clean=="MC")
 pt <- subset(leg_did, leg_did$party_clean=="PT")
 pvem <- subset(leg_did, leg_did$party_clean=="PVEM")
 pes <- subset(leg_did, leg_did$party_clean=="PES")
 panal <- subset(leg_did, leg_did$party_clean=="PANAL")
 other <- subset(leg_did, leg_did$party_clean=="Other")
 
 pan_m <- lm(sna_related ~ reelection, data=pan) 
 pri_m <- lm(sna_related ~ reelection, data=pri)
 mor_m <- lm(sna_related ~ reelection, data=mor)
 prd_m <- lm(sna_related ~ reelection, data=prd)
 mc_m <- lm(sna_related ~ reelection, data=mc)
 pt_m <- lm(sna_related ~ reelection, data=pt)
 pvem_m <- lm(sna_related ~ reelection, data=pvem)
 pes_m <- lm(sna_related ~ reelection, data=pes)
 panal_m <- lm(sna_related ~ reelection, data=panal)
 other_m <- lm(sna_related ~ reelection, data=other)
 
 # Put model estimates into temporary data.frames:
 pan_d <- data.frame(Variable = rownames(summary(pan_m)$coef),
                     Coefficient = summary(pan_m)$coef[, 1],
                     SE = summary(pan_m)$coef[, 2],
                     party = "PAN")
 pri_d  <- data.frame(Variable = rownames(summary(pri_m)$coef),
                      Coefficient = summary(pri_m)$coef[, 1],
                      SE = summary(pri_m)$coef[, 2],
                      party = "PRI")
 mor_d <- data.frame(Variable = rownames(summary(mor_m)$coef),
                     Coefficient = summary(mor_m)$coef[, 1],
                     SE = summary(mor_m)$coef[, 2],
                     party = "MORENA")
 prd_d  <- data.frame(Variable = rownames(summary(prd_m)$coef),
                      Coefficient = summary(prd_m)$coef[, 1],
                      SE = summary(prd_m)$coef[, 2],
                      party = "PRD")
 mc_d  <- data.frame(Variable = rownames(summary(mc_m)$coef),
                     Coefficient = summary(mc_m)$coef[, 1],
                     SE = summary(mc_m)$coef[, 2],
                     party = "MC")
 pt_d  <- data.frame(Variable = rownames(summary(pt_m)$coef),
                     Coefficient = summary(pt_m)$coef[, 1],
                     SE = summary(pt_m)$coef[, 2],
                     party = "PT")
 pvem_d  <- data.frame(Variable = rownames(summary(pvem_m)$coef),
                       Coefficient = summary(pvem_m)$coef[, 1],
                       SE = summary(pvem_m)$coef[, 2],
                       party = "PVEM")
 pes_d  <- data.frame(Variable = rownames(summary(pes_m)$coef),
                      Coefficient = summary(pes_m)$coef[, 1],
                      SE = summary(pes_m)$coef[, 2],
                      party = "PES")
 panal_d  <- data.frame(Variable = rownames(summary(panal_m)$coef),
                        Coefficient = summary(panal_m)$coef[, 1],
                        SE = summary(panal_m)$coef[, 2],
                        party = "PANAL")
 other_d  <- data.frame(Variable = rownames(summary(other_m)$coef),
                        Coefficient = summary(other_m)$coef[, 1],
                        SE = summary(other_m)$coef[, 2],
                        party = "Other")
 
 
 
 # Combine these data.frames
 coefs <- data.frame(rbind(pan_d, pri_d, mor_d, prd_d, mc_d, pt_d, pvem_d, pes_d, panal_d, other_d))
 
 ## Specify the width of your confidence intervals
 interval <- -qnorm((1-0.95)/2)  # 95% multiplier
 
 coefs$lower <- coefs$Coefficient - coefs$SE*interval
 coefs$upper <- coefs$Coefficient + coefs$SE*interval
 
 coefs <- subset(coefs, coefs$Variable=="reelection")
 
 # Order factors
 coefs$party <- factor(coefs$party, levels = c("PAN",
                                               "PRI", 
                                               "PRD", 
                                               "MORENA",
                                               "MC",
                                               "PT",
                                               "PVEM",
                                               "PANAL",
                                               "PES",
                                               "Other"))
 
 sna <- ggplot(data = coefs, aes(x = party, y = Coefficient)) + 
   geom_point() + 
   geom_errorbar(data = coefs, aes(ymin=lower, ymax=upper), width=0) + 
   theme_light() +
   xlab("") + 
   ylab("Estimate") + 
   ggtitle("SNA-related sponsorship") + 
   #ylim(-0.8, 1) + 
   theme_light() + 
   geom_hline(yintercept=0, colour = "grey60", linetype = 2) +
   theme(legend.position="bottom",
         legend.justification="center",
         panel.grid.major = element_blank(),
         panel.grid.minor = element_blank(),
         legend.title= element_blank(),
         axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
         plot.title = element_text(hjust = 0.5),
         plot.margin=unit(c(1, 1 , .8, 1.5),"cm"),
         axis.title.x.bottom = element_text(margin = margin(15, 0, 1, 0)),
         axis.title.y.left = element_text(margin = margin(0, 15, 0, 1)))
 
 
 
 ## ACIs
 pan_m <- lm(acl_total ~ reelection, data=pan) 
 pri_m <- lm(acl_total ~ reelection, data=pri)
 mor_m <- lm(acl_total ~ reelection, data=mor)
 prd_m <- lm(acl_total ~ reelection, data=prd)
 mc_m <- lm(acl_total ~ reelection, data=mc)
 pt_m <- lm(acl_total ~ reelection, data=pt)
 pvem_m <- lm(acl_total ~ reelection, data=pvem)
 pes_m <- lm(acl_total ~ reelection, data=pes)
 panal_m <- lm(acl_total ~ reelection, data=panal)
 other_m <- lm(acl_total ~ reelection, data=other)
 
 # Put model estimates into temporary data.frames:
 pan_d <- data.frame(Variable = rownames(summary(pan_m)$coef),
                     Coefficient = summary(pan_m)$coef[, 1],
                     SE = summary(pan_m)$coef[, 2],
                     party = "PAN")
 pri_d  <- data.frame(Variable = rownames(summary(pri_m)$coef),
                      Coefficient = summary(pri_m)$coef[, 1],
                      SE = summary(pri_m)$coef[, 2],
                      party = "PRI")
 mor_d <- data.frame(Variable = rownames(summary(mor_m)$coef),
                     Coefficient = summary(mor_m)$coef[, 1],
                     SE = summary(mor_m)$coef[, 2],
                     party = "MORENA")
 prd_d  <- data.frame(Variable = rownames(summary(prd_m)$coef),
                      Coefficient = summary(prd_m)$coef[, 1],
                      SE = summary(prd_m)$coef[, 2],
                      party = "PRD")
 mc_d  <- data.frame(Variable = rownames(summary(mc_m)$coef),
                     Coefficient = summary(mc_m)$coef[, 1],
                     SE = summary(mc_m)$coef[, 2],
                     party = "MC")
 pt_d  <- data.frame(Variable = rownames(summary(pt_m)$coef),
                     Coefficient = summary(pt_m)$coef[, 1],
                     SE = summary(pt_m)$coef[, 2],
                     party = "PT")
 pvem_d  <- data.frame(Variable = rownames(summary(pvem_m)$coef),
                       Coefficient = summary(pvem_m)$coef[, 1],
                       SE = summary(pvem_m)$coef[, 2],
                       party = "PVEM")
 pes_d  <- data.frame(Variable = rownames(summary(pes_m)$coef),
                      Coefficient = summary(pes_m)$coef[, 1],
                      SE = summary(pes_m)$coef[, 2],
                      party = "PES")
 panal_d  <- data.frame(Variable = rownames(summary(panal_m)$coef),
                        Coefficient = summary(panal_m)$coef[, 1],
                        SE = summary(panal_m)$coef[, 2],
                        party = "PANAL")
 other_d  <- data.frame(Variable = rownames(summary(other_m)$coef),
                        Coefficient = summary(other_m)$coef[, 1],
                        SE = summary(other_m)$coef[, 2],
                        party = "Other")
 
 # Combine these data.frames
 coefs <- data.frame(rbind(pan_d, pri_d, mor_d, prd_d, mc_d, pt_d, pvem_d, pes_d, panal_d, other_d))
 
 ## Specify the width of your confidence intervals
 interval <- -qnorm((1-0.95)/2)  # 95% multiplier
 
 coefs$lower <- coefs$Coefficient - coefs$SE*interval
 coefs$upper <- coefs$Coefficient + coefs$SE*interval
 
 coefs <- subset(coefs, coefs$Variable=="reelection")
 
 # Order factors
 coefs$party <- factor(coefs$party, levels = c("PAN",
                                               "PRI", 
                                               "PRD", 
                                               "MORENA",
                                               "MC",
                                               "PT",
                                               "PVEM",
                                               "PANAL",
                                               "PES",
                                               "Other"))
 
 aci <- ggplot(data = coefs, aes(x = party, y = Coefficient)) + 
   geom_point() + 
   geom_errorbar(data = coefs, aes(ymin=lower, ymax=upper), width=0) + 
   theme_light() +
   xlab("") + 
   ylab("Estimate") + 
   ggtitle("ACI sponsorship") + 
   #ylim(-0.8, 1) + 
   theme_light() + 
   geom_hline(yintercept=0, colour = "grey60", linetype = 2) +
   theme(legend.position="bottom",
         legend.justification="center",
         panel.grid.major = element_blank(),
         panel.grid.minor = element_blank(),
         legend.title= element_blank(),
         axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
         plot.title = element_text(hjust = 0.5),
         plot.margin=unit(c(1, 1 , .8, 1.5),"cm"),
         axis.title.x.bottom = element_text(margin = margin(15, 0, 1, 0)),
         axis.title.y.left = element_text(margin = margin(0, 15, 0, 1)))
 
 
 ggarrange(aci, sna, nrow=2, ncol=1)
 
 
 
 
 